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AN ASSESSMENT OF PREWHITENING IN ESTIMATING POWER SPECTRA 
OF ATMOSPHERIC TURBULENCE AT LONG WAVELENGTHS 

Samuel R. Keisler* and Richard H. Rhyne 
Langley Research Center 


SUMMARY 

An assessment is given of the effects of prewhitening on 
determination of power spectra of atmospheric turbulence at long 
wavelengths. In a computer experiment, a synthetic time history 
was generated by combining sine waves of random phase and fre- 
quency (but uniformly distributed over the frequency range of 
interest) with amplitudes adjusted to produce a power spectrum 
of known shape approximating that of atmospheric turbulence. The 
synthetic time history was then used to assess bias errors in 
power spectra computed by three different algorithms implemented 
by the fast Fourier transform. 

Prewhitening is unnecessary when using the narrow "spectral 
windows" required for determining power spectra of atmospheric 
turbulence below the "knee" frequency, or at very long wavelengths. 
Several wide spectral window oases are also included to assess the 
effect of first-difference prewhitening on data where the first 
spectral estimate is above the knee frequency. 

INTRODUCTION 

The generally accepted mathematical model of the power spec- 

A ^ 

density function for atmospheric turbulence (the Von Karman 

*SDC Integrated Services, Inc. 
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model; see ref. 1) describes the amplitude of the spectrum as 
nearly constant at very low frequencies, varying inversely with 
the five-thirds power of frequency at the higher frequencies. 

The transition between these frequency regions is often referred 
to as the "knee” and occurs at a frequency that is dependent on 
the integral scale of turbulence (L). The knee frequency, and 
therefore the scale of turbulence, is significant with respect to 
the calculated motion and load responses of aircraft (ref. 1). 

The adequacy of the Von Karman spectrum as a model of atmo- 
spheric turbulence has not been assessed sufficiently. Although 
the shape of the spectrum at high frequencies has been shown to 
agree with measured turbulence spectra, there is still uncertainty 
about the representation by the Von Karman spectrum at frequencies 
associated with and below the knee of the theoretical spectrum. 
This is caused, to a large extent, by difficulties in accurately 
determining values in the low frequency region of spectra of atmo- 
spheric turbulence from flight measurements. These difficulties 
include (1) obtaining adequate length samples, (2) Inaccuracies 
in measurements of airplane motions that are required to extract 
the turbulence velocities, and (3) errors introduced by the con- 
version of the measured time histories of the velocities to power 
spectral density functions. Efforts are underway to reduce errors 
from the first two sources. This report is concerned with the 
third problem cited. 

A significant source of error in the determination of experi- 
mental spectra is the bias error associated with data processing. 
The bias error is negligible for spectra of essentially constant 
amplitude, but it can produce an unacceptably large distortion in 
spectra exhibiting large changes in amplitude. Therefore, the 
high-frequency region of the turbulence spectrum, where the ampli- 
tude varies inversely with the five-thirds power of frequency, is 
susceptible to bias error. Press and Tukey (ref. 2) introduced a 
conditioning process called prewhitening to reduce the distortion 
of turbulence spectra substantially (that distortion caused by 
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data processing) at the higher frequencies since the higher fre- 
quency range was the region of particular interest at the time. 

The prewhitening procedure, which has become routine, renders the 
high-frequency range of the spectrum to nearly constant amplitude 
during processing. The effect is then removed (postdarkened ) for 
final data presentation. If, however, the actual turbulence spec- 
trum in the low-frequency region is nearly constant, as suggested 
by the Von Karman spectrum model, the prewhitening process results 
in a bias distortion of the spectrum in this low-frequency range. 
Current interest in Improved definition of measured spectra at low 
frequencies prompted a reassessment of the use of the prewhitening 
process and other aspects of the bias error problem. 

The present computer study^ is similar in approach to the 
study described in reference 3 in that spectra are estimated by 
the digital processing of a synthetic random time function that 
simulates the velocities of atmospheric turbulence. In contrast 
to reference 3> which examined the effects of various methods of 
prewhitening, the present study compares spectra obtained with and 
without prewhitening, with emphasis on the effects at low frequen- 
cies. The method used in this study to generate the synthetic ran- 
dom time function differs from the method used in reference 3; the 
reasons for this different method are discussed in a later section. 
The data processing methods consist of three commonly used fast 
Fourier transform (FFT) algorithnjs: (1) Blackman-Tukey , 

(2) ensemble averaging, and (3) frequency averaging. 

SYMBOLS 

Values are given in both SI and U.S. Customary Units. The 
measurements were made in U.S. Customary Units. 

^Recommended by Mr. Allan Piersol (Bolt, Beranek, and Newman, 
Inc., Canoga Park, CA). 
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B, 


equivalent spectral bandwidth, Hz 


B 


®k' 


dimensionless equivalent spectral bandwidth, 
2 ttL 



f frequency, Hz 

smoothed power spectral density estimate at frequency 

i Af 


G» 


raw power spectral density estimate at frequency 
i Af 


H(f ) 


transfer function 


i,j,k,n index, integer variable 


K 


number of points, including zeros, transformed 
by fast Fourier transform (power of 2) 


dimensionless wave number. 


2TrfL 

V 


L integral scale of turbulence 

* number of nonnegative lags of time-sampled autocorrelation 

function 

M number of contiguous raw power spectral density estimates 

which are used to compute each point of frequency aver- 
aged power spectral density; or. number of signal seg- 
ments used in computing ensemble average power spectral 
;denslty 
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integer constant, i.e., final integer in function 
sequence 

time length of function (or segment of function) 
being transformed 

frequency-sampled spectral estimate at frequency 

i Af 

turbulence sample length, sec 

time, sec 

velocity 

random time function 

time function sampled at time i At 

prewhitened random time function 

prewhitened time function sampled at time i At 

frequency interval, Hz 

time sampling interval, sec 

phase angle, rad 

lag time 

maximum lag time 


power spectrum as function of f 



# power spectrum as function of k’ 

true Dryden spectrum undistorted by bias error 
power spectrum of postdarkened time function 
*y power spectrum of prewhitened time function 

0 ) angular or circular frequency, 2irf 

0 ) ^ circular frequency Integration variable 

APPROACH 

Since the primary purpose of this study is to assess the value 
of prewhitening in reducing bias error in estimated spectra at low 
frequencies, a review of bias error is appropriate. This review 
includes: (1) properties related to the general digital processing 

of a random time function, (2) effects at high and low frequencies, 
and (3) a brief description and discussion of the prewhitening pro- 
cess. The review is followed by descriptions of the three digital 
algorithms for estimating power spectral densities and their rela- 
tive bandwidth properties. The properties of a synthetic random 
time function needed to represent atmospheric turbulence velocities 
are then discussed, and the two functions used in numerical calcu- 
lations are described. The conditions selected for estimating 
power spectral densities are stated and results are discussed. 

BIAS ERROR 

Each of the power spectral density algorithms considered here 
uses the Fourier transform of a time function of finite length. The 
following discussion is equally applicable to the Fourier transform 
of a finite length autocorrelation function as well as a finite 
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length time function. Further details are given on page IVC;19 of 
reference 2. The finite time length of a time function affects the 
measured spectral function and can be considered the product of an 
infinite length function (given the subscript "true”) and a "boxcar” 
function (fig. 1(a)). The measured and the true time functions are 

-P P 

equal when — < t < -. The measured time function is equal to zero 
2 2 

for all other values of time. The Fourier transforms of the two 
time functions are related by the convolution integral 


T, ( <»> ) = 

meas ' 


♦true(“l^P 



p 

sin 

(o) - U 1 ) — 

2 


do) 


(o) - 0) . )- 


( 1 ) 


where *true^“^ "true" spectrum for limit p ->• ® and the 

sin X 

— — term is the spectral "window." (See fig. 1(b).) The differ- 
ence between *meas^“^ *true^“^ defined as the bias error. 

(This example is applied to the boxcar window; however, the same 
relation would hold for other windows.) Convolution can be under- 
stood heuristically by considering its effect on an arbitrary spec- 
tral estimate *rneas^“^> according to the following procedure: 

Shift the spectral window so that it is centered at u > ; form -the 
product of the true spectrum *true shifted window; then 

*meas^“^ is the integral of this product. 
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The strongest objection to the boxcar window concerns "leakage 
through the side lobes" (fig. 1(b)); the side lobes are those por- 


1 


2 

2 

3 

tions of the window between - 

and 

-, between 

— and 

“, and 

P 


P 

P 

P 

1 

-1 

1 


the "main 

so forth. [The portion between 
\ 

— 

and “ is 

called 

P 

P 




lobe.’y When the process of convolution is kept in mind, it can be 
seen that any portion of the true spectrum which is multiplied by 
a side-lobe peak during convolution can have a significant effect 
upon the estimate. 

This deficiency in the boxcar window is largely overcome in 
windows, such as the Hann window (ref. 4), shown in figures 2(a) 
and 2(b), which taper the ends of the function. The Hann side 
lobes are greatly reduced in relation to the boxcar side lobes. 

Figure 3 Illustrates a typical bias error encountered in esti- 
mating the power spectrum. The true spectrum together with the 
shifted spectral window are shown with linear scale on each axis. 

The true spectrum *true^“^ shown in figure 3 is an approximation 
of the Von Karman turbulence model, and is discussed later in more 
detail. It can be visualized that the magnitude and sign of the 
bias error depend primarily upon the broadness and general shape of 
the main lobe of the spectral windown and its position along the 
frequency axis with respect to the true spectrum. In general, the 
peak of the true spectrum is underestimated, and all other values 
overestimated . 

Each of the three power spectral density algorithms is asso- 
ciated with a different variant of the fundamental spectral window 
described above. The differences among these variants are affected 
very little by the time window chosen, e.g., Hann, Hamming, or Parzen 
(see ref. 4). For this reason, only the Hann time window is consid- 
ered here. 
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PREWHITENING 


The prewhitening process consists of filtering the random time 
function so as to produce a nearly flat spectrum at frequencies 
above the knee. (The purpose, as stated earlier, is to reduce or 
eliminate the bias error caused by the window convolution illus- 
trated in fig. 3.) The spectral estimates are then ”postdarkened” 
to correct for the initial prewhitening. The most commonly used 
prewhitening process is the first-difference method 

y(t ) = x(t ) - x(t - At ) 
or in digitized sample form 

Yi = Xf - x^_^ (2) 

The spectrum estimated from the prewhitened function by a particu- 
lar digital algorithm is then postdarkened as follows before final 
data presentation 


^x(f) 



*y(f) 


(3) 


where 


1 1 
|H(f)|2' 2(1 — cos 2irf At) 

(See ref. 2, or appendix E of ref. 8.) 

The true spectrum »^(f) of figure 3 is shown in figure 4 
as #y(f) after the first-difference prewhitening filter has 
been applied, i.e., |H(f)|^. The knee frequency 

of figure 3 is indicated by an arrow on the abscissa of figure 4. 
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(The frequency scale in fig. 4(b) is exaggerated in comparison 
to the scale of fig. 4(a) to show more clearly the low-frequency 
region.) It is apparent that the bias error will be posi.tive at 
the knee frequency and below, and will become larger as the spec- 
tral window becomes broader. At frequencies well above the knee, 
as the prewhitened spectrum becomes more nearly flat, the bias 
error will become negligible, even though the spectral window is 
quite broad. 


POWER SPECTRAL DENSITY ALGORITHMS 

The turbulence time histories are sampled in time at interval 
At. The time function is denoted by and is the value of the 

function sampled at time i At (where i = 1, 2, . . ., N); then 

discrete spectral estimates s^' corresponding to harmonic frequen- 

• ^ 

cies 1 Af are computed by the EFT where i = 0, 1,2, . . ., - 

\ 2 

1 \ 

and Af = . The time length of the function (or segment of 

N At/ 

function) being transformed is P = N At; so 

1 1 

Af = - = (4) 

P N At 


where Af is the frequency spacing of the original "raw” 
estimates. Each of the power spectral density algorithms 
is now examined together with its associated spectral window. 

Blackman -Tukey 

The technique for computing the autocorrelation function 
with Fourier transforms rather than with lag products as used 
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in the past is given in reference 4 (pp. 165 and 166). The 
Blackman-Tukey FFT power spectral density algorithm is outlined 
briefly in the following eight steps. 

(1) Let K, which must be a power of two, be the number of 
points being transformed. Let *, be the number of nonnegative 
lags required for the autocorrelation function. Only K - 

data points can be used, since the last i values of the points 
being transformed, at least, must be zeroed to avoid the circular 
autocorrelation effect (ref. 4, pp. 123 and 124). The mean of 
the data function is removed and the correct number of zeros 
is added. 

(2) Optional prewhitening: the first difference filter is 

applied (eq . (2) ) . 

(3) The FFT is computed to yield the complex amplitude 
spectrum. 

(4) The ’’raw" power spectrum is computed as the scaled 
square modulus of the result of step (3)* 

(5) The inverse FFT is computetl to produce the autocorrela- 
tion function. 

(6) The Hann window is applied. 

(7) The FFT is computed to obtain the power spectrum, which 
is then scaled to obtain the power spectral density. 

(8) Postdarkening (if option (2) is selected) is applied as 
indicated by equation (3). 

The autocorrelation function computed in step (5) is defined 

“"^m = = "^m where t is lag time and "^m = “ '')^t. 

The quantity a is adjusted to truncate the autocorrelation func- 
tion at an appropriate maximum lag Tj^ which results in the desired 
spectral window width B^. For convenience in the FFT procedure, 

«, - 1 should also be a power of two. If is set equal to 

P 1 

— in equation (4), Af □ — — . When the Hann window shown in -fig- 


ure 2 is used in step (6), then the spectral windows for two adjacent 



spectral estimates at frequencies 1 Af and (i + 1)Af will 
be as shown in figure 5. The mechanics of the Blackman-Tukey 
algorithm are such that the frequency spacing of the final power 

1 

estimates Af is always , and the effective spectral window 


1 

width B is — or 2 Af. A rather large degree of window 
® Tm 

overlap for adjacent spectral estimates is shown in figure 5. 

It is apparent, therefore, that such estimates are not strictly 
independent. For this reason, it is sometimes advocated that 
every other estimate be discarded when employing this algorithm. 

Consider now the estimate Sq at zero frequency. The zero 
frequency component of the signal (i.e., mean) was removed in 
step ( 1 ) . Therefore it appears that Sq is an estimate of the 
power at positive frequencies. For this reason, Sq is (somewhat 

Af 

arbitrarily) displayed at frequency — in this paper. 


Ensemble Averaging 

This algorithm is described in detail in reference 5. A 
brief outline follows. 

(1) The time function is partitioned into M segments of 
equal length, with the number of points in each segment a power 
of two. Steps (2) to (6) must be performed independently for 
each segment. 

(2) The mean of the function is removed. 

(3) (Optional prewhitening) The first difference filter is 
applied . 

(4) The side-lobe suppression window is applied. 

(5) The FFT is computed to obtain the complex amplitude 
spectrum. 
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(6) The power spectral density estimates are computed as the 
scaled square modulus of the amplitude spectrum. 

(7) The M segment power spectral density functions are 
averaged . 

(8) (Optional) Postdarken. 

In reference 5, the observation is made that the expected 
value of the spectrum at each point is the convolution of the true 
spectrum with the square of the spectral window if the time window 
is even (a proof of this is given in ref. 6). All commonly used 
windows (including the Hann window) are even. 

The squared Hann window is shown in figure 6. It is the 
effective window in the limit as the length of the signal, and 
therefore the number of segments used, increases without bound. 

For this reason, the squared spectral window will be referred to 
here as the "expected window." 

Figure 6 shows that the ensemble average window bandwidth 
Bg, that is, the width of the main lobe at the half-power point. 


is approximately 1 
segment; then when 
power estimates is 


5 Af. Let - 
M 


- = N ( from 
M 


Af, equal to 


be the number of points per 


eq. (4)), the spacing of the 


M 

or 

K At 


B 

e 


1.5 Af 


1 .5M 

K At 


Note that in this figure contiguous windows overlap less 
than for the Blackman-Tukey algorithm; therefore, contiguous 
ensemble average estimates are statistically more independent 
of each other. 
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Frequency Averaging 


This algorithm is a very simple modification of the ensemble 
averaging algorithm. 

(1) Perform steps (2) to (6) of the ensemble averaging method 
for the entire time function consisting of K points. Denote the 

K 

resulting spectrum where i = 0, 1, . . • » 

(2) The smoothed power spectral density estimates, here 
denoted Gj, are computed as the means of M contiguous values of 

as follows: 




G* +G* _+ . . . +G* . 

k+1 k+2 k+M 


M 


and 


j 



(j - DM 


(3) (Optional) Postdarken as above. 

The expected window is now the sum of M contiguous windows, 
each as shown in figure 6 (see ref. 7, p. 296). The frequency 
average expected window for M = 16 is shown in figure 7; the 
estimate is displayed at the midpoint of the window, i.e., the jth 

(2j - DM + 1 

point of the spectrum is displayed at frequency . 

2K At 

The final frequency spacing of the power estimates Af is then 

M M + 0.5 

, and the bandwidth B is approximately . As can be 

K At ® K At 


n 



seen in figure 7, the expected window is now essentially rectan- 
gular, with virtually no overlap of contiguous windows. 

In comparison to ensemble averaging, frequency averaging has 
this additional advantage: for the same frequency spacing of the 

final power estimates (and approximately the same bandwidth 
the lowest frequency analyzed is lower by approximately a factor 


M + 1 

of two. The lowest frequency estimate appears at for fre- 


quency averaging, as opposed to 


M 

K At 


2K At 

for ensemble averaging. 


SYNTHETIC RANDOM TIME FUNCTIONS 


For this specific study, the synthetic random time function 
must possess the significant characteristics of atmospheric tur- 
bulence velocities. These include randomness, a nearly Gaussian 
probability density function, and a power spectral density repre- 
sentative of that for turbulence. The power spectral density was 
chosen to be the Dryden spectrum (see ref. 8) to parallel the study 
of reference 3. The differences between the Dryden and Von Karraan 
spectra are not considered to be significant for the present pur- 
poses. The Dryden spectrum is described by 


i (k’ ) 


1 1 + 3k’ 

ir 


(1 + k'^) 


2-2 


(5) 


where 


2irfL 

k’ = 

V 
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(sometimes referred to as "dimensionless wave number"), and where 


L integral scale of turbulence 

V velocity of aircraft 

f frequency, Hz 

In terms of temporal frequency 


$(L,V,f) 



( 6 ) 


The synthetic function, of practical necessity, is synthesized 
from a finite set of numbers and therefore, cannot completely 
describe a power spectral density that is continuous in frequency. 
Consequently, it is necessary to analyze any particular synthetic 
signal for adequacy in this respect. This was done for the signal 
used in reference 3 (designated herein as Signal B) and details are 
given in the appendix. Signal B is a model of a process whose fre- 

i 

quency components are harmonics of the fundamental frequency -. 

It is therefore not a good representation of atmospheric turbulence, 
whose frequency components are distributed continuously over all 
frequencies. The frequency components of Signal B coincide pre- 
cisely with the nulls of the boxcar window and with the center of 
the main lobe. Consequently, Signal B can produce no bias error 
with the boxcar window. This effect is considered unrealistic with 
respect to simulating atmospheric turbulence. 
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To model the continuous frequency distribution of turbulence 
more accurately, an alternate synthetic time function, designated 
Signal A, was derived for use in this study. Details concerning 
this signal are also given in the appendix. 

Signal A = ^ H(f) sin ( 2 irfjjn At + e^) 

where both fp and 0 j^ are randomly picked values in each of 
1250 frequency intervals. Thus, spectral values do not neces- 
sarily coincide with boxcar window nulls, and the approximation 
of atmospheric turbulence is improved. A short segment of the 
time function is shown in figure 8 ; its appearance is very much 
like that of actual turbulence time-history measurements. The 
cumulative distribution function, shown in figure 9 , is very 
nearly Gaussian, as desired, and as expected from the procedure 
used. 


CONDITIONS CONSIDERED FOR SPECTRAL ESTIMATION 

To determine a practical range of bandwidth Bg which 
should be investigated for bias error effects, constraints 
imposed by sample-length limitations of actual turbulence data 
must be considered. It would, of course, be desirable to have 
the spectral estimates very closely spaced (small Bg) in the 
frequency region approaching zero (see fig. 3 ) to define a knee, 
if one exists. As the integral scale of turbulence L becomes 
larger in the mathematical model, the knee frequency moves closer 
to zero frequency, and in fact might be very close to zero fre- 
quency for certain meteorological conditions. For example, if L 
is 762 m (2500 ft) in the Von Karman turbulence model (a present 
turbulence design requirement; see ref. 9 ) and the velocity of 
the airplane is I83 m/s (600 ft/sec), the knee peak appears at 
0.017 Hz. Larger L values result in lower knee frequencies. 

As the bandwidth is made smaller (for a given sample length 
of actual turbulence data), the random fluctuations of the power 

17 
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estimates become larger, eventually making it impossible to 
determine whether a knee is present. Past experience and a cal- 
culation of the magnitude of the random power fluctuations based 
upon so-called "statistical degrees of freedom" (see ref. 7) indi- 
cate that approximately 30 degrees of freedom is as low as can be 
tolerated . 

According to reference 7 (p. 219), 

Degrees of freedom (d.o.f.) = 2B T (7) 

where is bandwidth in Hz, and T is total sample length in 

seconds. Past experience with turbulence data collection indi- 
cates that usable sample lengths of more than 10 minutes duration 
are not readily obtainable, particularly at the higher altitudes. 
If we assume the maximum practical sample length to be 10 minutes 
(or 600 sec), and that 30 d.o.f. can be tolerated, then B^ is 
0.025 Hz. 

Based on the preceding discussion, indications are that 
although narrower bandwidths for better spectral resolution at 
very low frequencies would be desirable, bandwidths smaller than 
about 0.025 Hz are not practical because of sample length limita- 
tions. The present bias error investigation will therefore be 
confined to bandwidths on the order of 0.025 Hz and greater. (It 
should be noted that random error is only of indirect concern in 
this study, and that it has, in fact, been minimized by the syn- 
thetic time functions so as not to obscure the bias error under 
study. ) 

Since this study was initiated by following the procedure of 
reference 3> the numbers actually used in equation (6) to generate 
the synthetic signals were L = 191 m (627 ft) and V = I 83 m/s 
(600 ft/sec) to conform with reference 3. The synthetic time 
function values were spaced at an Interval of At = 0.04 sec, and 
this spacing provided a Nyqulst frequency of the resulting power 
spectral density of 12.5 Hz, also in accordance with reference 3. 
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A more general, and possibly more useful, procedure would 
have been to generate the spectral results in terms of dimension- 
less wave number k' (see eq. (5)) rather than in terms of tem- 
poral frequency. The correct shape of the spectrum is all that is 
required, however, and the exact location of the knee peak along 
the frequency axis will not affect, the magnitude of the bias error 
obtained. If the location of the knee peak is known in advance, 
and/or if an estimate of bias error for a specific L, B^, or V 
at a particular frequency relative to the knee peak is needed, then 
k' units can be helpful. A k' scale has therefore been included 
on all the spectral plots presented. Also B is given in k* 
units as well as in Hz, and is then designated B 

^ k ^ 

The synthetic time function was processed as if it were mea- 
sured data. Spectra were estimated by the three algorithms with 
and without first-difference prewhitening. As indicated, the min- 
imum bandwidth desired was approximately 0.025 Hz. The nearest 
bandwidths readily obtainable with the three algorithms were 
Bg = 0.0244, 0.0366, and 0.0260 Hz for Blackman-Tukey , ensemble 

averaging, and frequency averaging for B. = 0.160, 0.240, and 

\ \ ®k ' 

0.171 ), respectively. These values were determined as discussed 

under power spectral density algorithms, by use of 1024 nonnega- 
tive lags for Blackman-Tukey, 15 time segments each containing 
1024 data points for ensemble averaging, and frequency averaging 
over 16 contiguous raw estimates obtained by transforming 15 871 
data points plus 513 zeros to make an even power of two. The 
Af frequency spacing was thus 0.0122 Hz for Blackman-Tukey and 
0.0244 Hz for ensemble averaging and frequency averaging. 

Additional comparisons between prewhitened and non-prewhitened 
spectra were made with the Blackman-Tukey algorithm for an inter- 
mediate Bg of 0.195 Hz ^Bg ^ = 1.28^ employing 128 lags, and for 
a very wideband Bg of 0.391 Hz or Bg ^ of 2.56 (64 positive 
lags). The wideband spectra simulated results obtained in earlier 
years before FFT was available. 
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RESULTS AND DISCUSSION 


In the following figures, the power spectral density is plot- 
ted along the vertical axis, and the frequency in Hz and in k’ 
units is plotted along the horizontal axis, both on logarithmic 
scales. The true power spectral density ♦.j.(f) is modified by 
a multiplication factor selected so that the area under *^(f) 
will be equal to the variance of the synthetic random time func- 
tion. The obtained *-p(f) is superimposed on the plot of each 
computed power spectral density. In each case, the "true” spec- 
trum is labeled "reference" and the spectrum calculated from the 
time history is labeled "computed.” 

The results for the minimum bandwidth case, presented in fig- 
ures 10, 11, and 12, indicate that all three algorithms define the 
knee adequately. There is an indication of bias error (a slight 
overestimation in the estimates nearest zero frequency) only in the 
spectrum obtained with the Blackman-Tukey algorithm (fig. 10). 

With all other parameters held fixed, first-difference pre- 
whitening and postdarkening were applied. In each c.ase the effect 
(see figs. 13, and 15) was an overestimation at the lowest fre- 

quencies and no significant improvement elsewhere. The curiously 
large difference between the amount of distortion obtained by using 
the Blackman-Tukey algorithm and the amount of distortion obtained 
by using the direct-transform algorithms needs an explanation. To 
rule out the possibility of error in the Blackman-Tukey portion of 
the computer program, the synthetic signal was processed through a 
lag-product Blackman-Tukey program (see ref. 10), and the results 
were in close agreement with those obtained from using the FFT 
version . 

The first-difference prewhitened true Dryden spectrum is shown 
in figure 16; it is the product of the first-difference transfer 
function squared |H(f)|^ and the true Dryden spectrum of equa- 
tion (6). Since this curve is concave upward in the frequency range 
of the first several points of the spectrum, the estimates at these 
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points should contain positive bias errors arising from the window 
convolution described earlier. The postdarkener function (eq. (3)) 
in this frequency range is quite large, going to infinity at zero 
frequency; therefore, the postdarkening operation could transform 
the positive bias errors mentioned earlier into a curve similar to 
figure 13 or to figures 14 and 15. 

To determine the precise amount of bias error caused by window 
convolution alone, a digital simulation of the Blackman-Tukey win- 
dow convolution was performed as follows: The Hann spectral window, 

1 

where Af = — is the same as the Af for figure 13 (see fig. 2(b)), 
P 

was represented on the interval -8 Af to 8 Af as a 16 001-point 

t' 

sampled-data function. The prewhitened Dryden spectrum, including 
the negative-frequency plane, was also represented as a sampled- 

Af 

data function, with the same frequency spacing . The con- 

1000 

volution integral was then computed at frequencies i Af, where 
i = 1, 2, . . 10 from the vector dot product of the shifted 

window and the theoretical curve. These 10 values were then post- 
darkened. The maximum bias occurred at i = 1, the spectrum at 
this point being approximately 1.5 times the correct value. There- 
fore, window convolution bias alone is not sufficient to account 
for the Blackman-Tukey results of figure 13, although it probably 
accounts for a substantial portion of the error in the two direct- 
transform algorithm results (figs. 14 and 15). 

Two additional possible sources of error are now discussed. 

The autocorrelation function (see ref. 8, p. 13) indicates a high 
correlation between successive values of the time function (that 
is , at T = At ) . This implies that the application of the first- 
difference filter (eq. (2)), where the subtraction of successive 
points occurs, results in a loss of precision of the filter output 
relative to the precision of the input. (Points of nearly equal 
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magnitude are being subtracted.) Because such round-off error is 
random, it adds a certain amount of white noise to the prewhitened 
spectrum and thus contributes to the positive bias near zero fre- 
quency. However, figures 14 and 15 suggest that this contribution 
is evidently small. Such an effect would be equal in all three 
algorithms . 

The final source of error considered here is the computer round 
off error which occurs during calculation of the power spectral den- 
sity from the time function. It should be noted that the Blackman- 
Tukey algorithm requires more than twice the number of computations 
as the other two algorithms (two full-length transforms plus one 
short transform for Blackman-Tukey , as compared with one full-length 
transform for the other two algorithms). The lag-product Blackman- 
Tukey algorithm requires still more calculations than the FFT ver- 
sion does. This computer round-off error is also random and could 
contribute a flat error spectrum to the positive bias near zero fre- 
quency. Therefore, if the white noise spectrum for the Blackman- 
Tukey algorithm is greater than those for the other two algorithms, 
and is as large as the approximation in figure 16, it could con- 
tribute to window convolution bias error and could account for the 
effect observed in figure 13. This is probably the source of the 
large difference observed in figure 13 in comparison with the dif- 
ferences shown in figures 14 and 15. 

The results obtained when a very wideband window is employed 
with the Blackman-Tukey algorithm ^B^ = 0.391 Hz or ^ = 2.56^ 

are now examined. Windows of about this size, together with first- 
difference prewhitening, were generally used for processing turbu- 
lence data before the FFT was available. (The data presented here, 
although processed by the use of FFT, are mathematically equivalent 
to the previously used lag-product Blackman-Tukey method, as indi- 
cated earlier.) The spectra with and without prewhitening are shown 
in figures 17(a) and 17(b), respectively. The corresponding pair 
for an intermediate bandwidth of 0.195 Hz = 1.28^ are shown 

in figure 18. In the past, the lowest frequency estimate from the 
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Blackman-Tukey lag-product algorithm was assumed to be at exactly 
zero frequency and was discarded because the postdarkener filter 
function (eq. (3)) goes to infinity at zero frequency, and also 
because the turbulence measurements were known to be contaminated 
by instrumentation drift problems at frequencies close to zero. 

If this lowest frequency point (which is postdarkened and dis- 

Af 

played at — on these plots, as indicated earlier) is discarded, 

4 

the lowest frequency estimate then appears at 0.195 Hz in fig- 
ure 17(a) and agrees very well with the reference curve. The next 
power estimate, which appears at 0.39 Hz, slightly underestimates 
the reference curve. For the non-prewhitened result (fig. 17(b)), 
the corresponding points underestimate and overestimate the true 
spectrum to a somewhat greater extent, so that for this large band- 
width the prewhitened results are indeed better (when the "zero 
frequency” estimate is ignored). On the other hand, figure 18 
shows that when the bandwidth is decreased to 0.195 Hz, the non- 
prewhitened results are better, whether the lowest frequency power 
estimate is ignored or not. 

These intermediate bandwidth results indicate that in the past, 
the use of first-difference prewhitening could have obscured a pos- 
sible knee in the experimental data. (See fig. 18(a).) It should 
be recalled, however, that interest at the time was directed toward 
defining the higher frequency end of the spectrum (using wider band- 
widths) with relatively short turbulence sample lengths. 

As an illustration of the minimum sample length required for 
this bandwidth, assume that use of 30 d.o.f. in equation (7) does 
not result in intolerably large random power fluctuations. A B 

* 

of 1.28 produces a lowest power estimate at about the knee peak 
location, and requires a sample length of only about 77 sec for the 

L 

ratio of — assumed here. 

V 
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Trade-off factors for real data are such that analyzing power 
spectra to longer wavelengths (and thus lower frequencies) requires 
a smaller ^also ®e{^t)> which in turn requires longer sample 

lengths to suppress the random statistical fluctuations of the power 

estimates. As B becomes smaller than about 1.2, which it must 

w T 

to analyze longer wavelengths or adequately resolve a possible knee, 
first-difference prewhitening is no longer desirable and, in fact, 
increases the bias error if the true spectrum is similar in shape to 
the Von Karman or Dry den turbulence model. 

The practical procedure should therefore consist of choosing a 
minimum B^ based on the sample length available. A determination 
or "rule of thumb" in terms of B„ will not suffice, since it 

t 

must be assumed that L, the integral scale of turbulence, is not 
known at this point, nor can it be assumed that the Von Karman or 
Dryden turbulence model is applicable, for that matter. Obtaining 
sufficiently long data samples to suppress the random error, which 
is considerably greater than bias error at this point, poses a prob- 
lem. It has been suggested that the ensemble averaging algorithm 
could be used to alleviate the difficulty by combining segments 
of data obtained from repeated passes through the same turbulent 
area, or from different turbulent areas known to have been generated 
by the same meteorological process. Such a procedure might be ques- 
tionable, however, since it would be based on the assumption that 
the spectral characteristics of all the segments would be Identical. 

If it is necessary to use an intermediate size bandwidth (that 
is, around 0.1 to 0.2 Hz) because of sample length limitations, some 
special type of prewhitening other than first difference could be 
beneficial. The spectral shape and frequency location of the knee 
must then be estimated or known in advance. For example, long sam- 
ples in a specific meteorological condition could have previously 
been processed and would lend confidence to a particular shape and 
knee location. Care must be exercised, however, since a "wrong 
guess" would lead to greater distortion of the final result, rather 
than less. 
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If wide bandwidth windowing is necessary (that is, of 

* 

about 2.6 and larger) because of limitations in sample length, then 
first-difference prewhitening is required for accurate spectral 
estimates to be obtained, as shown by the results given in figure 17. 
It is also recommended that the "zero frequency” estimate obtained 
with the Blackman-Tukey algorithm be discarded, as was done in the 
past, especially in view of the large distortion obtained near zero 
with this algorithm when using first-difference prewhitening. The 
spectral knee could not, of course, be resolved in this case since 
the lowest frequency spectral estimate obtained would be above the 
knee location. 


CONCLUDING REMARKS 


A study was made of the effects of prewhitening on determina- 
tion of power spectra of atmospheric turbulence at long wavelengths. 

A synthetic time history was generated by combining sine waves of 
random phase and frequency, amplitudes being adjusted to produce a 
power spectrum of known shape approximating that of atmospheric tur- 
bulence. The synthetic time history was then used in a computational 
experiment to assess bias errors in power spectra computed with and 
without prewhitening. 

Results of this experiment show that for minimum bandwidths 
deemed practical for processing 10-minute data samples (equivalent 
spectral bandwidth of 0.025 Hz), the finite bandwidth bias errors, 
for each of three power spectral density algorithms implemented by 
the fast Fourier transform, are negligible and are a great deal 
smaller than the random type errors expected. Prewhitening is 
therefore not recommended when power spectral estimates are obtained 
by employing these narrow bandwidths. First-difference prewhitening 
in particular was shown to have an undesirable effect upon power 
estimates at the lowest frequencies, especially when used with the 
Blackman-Tukey algorithm. 
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Fir^st-dif f erence prewhitening is not recommended when inter- 
mediate size bandwidths (that is, around 0.1 to 0.2 Hz) are used, 
although some special type of prewhitening could be beneficial if 
the spectrum shape and knee location were adequately known. First- 
difference prewhitening is, however, recommended for power estimates 
above the ’’knee” frequency for bandwidths of 0.3 Hz and greater. 

For model studies of spectral bias errors of a random time 
series, the function-generation technique used here is favored over 
the technique used by Otnes, Nathans, and Enochson ( AFFDL-TR-69-1 1 ) . 
The random spectral errors observed can be made arbitrarily small 
for a given function length by summing a sufficiently large number 
of random-frequency, random-phase sinusoids. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
August 11, 1976 
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APPENDIX 


SYNTHETIC RANDOM TIME FUNCTIONS 

Synthetic random time functions are often used in lieu of 
random time functions representing physical random processes. 

For example, aircraft flight simulators include responses to 
atmospheric turbulence velocities that are approximated by ran- 
dom time functions generated by suitably filtered electronic 
noise signals or by digital computer programs. This approach has 
the advantage of providing (a) arbitrarily chosen intensities, 

(b) stationarity if desired, and (c) simulation of the effect of 
various integral scales of turbulence. 

In the construction of digitally generated random time func- 
tions, care must be taken to assure that the significant features 
of the physical random process are modeled adequately. The digi- 
tal function, of course, must be constructed from a finite number 
of randomly chosen elements. There is the possibility that a par- 
ticular choice of elements may not satisfactorily model the physical 
process for certain purposes even though the probability distribu- 
tion and spectral shape are correct and the random time function 
appears to be representative of the physical process. This appen- 
dix presents an example of such an approximation (designated Sig- 
nal B) for atmospheric turbulence with respect to the assessment 
of effects of various spectral windows. An alternative approxima- 
tion, described in the text of this report (designated Signal A) 
is also discussed, with the generation technique described in more 
detail . 


Signal B 

Signal B was generated as follows, in accordance with 
reference 3: 
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(1) 16 384 random numbers were generated. These numbers were 
independent, Gaussian, with zero mean and unit variance. The At 
was selected to be 0.04 sec. 

(2) The complex FFT of these data was multiplied by the desired 
transfer function H(f), the square root of equation (6), with 

L = 191 m (627 ft) and V = 183 m/s (600 ft/sec). 

(3) The inverse FFT was then taken. 

The resulting function from step (2) is defined only at har- 
monically related frequencies n Af where n=0, 1, 2, . . ., 

8192, and where Af = ^ . The inverse transform of each of 

16 384 At 

these is a pure harmonic sinusoid with frequency n Af; therefore 
Signal B is the sum of the harmonic sinusoids. The spectrum of 
the model differs significantly from that for atmospheric turbu- 
lence since the turbulence spectrum is continuous and contains 
nonharmonically related components. 

The analysis of Signal B by the boxcar spectral window 
(fig. 1(b)) in equation (1) yields no bias error at all as a 
result of the coincidence of nulls in the spectral window with 

n 

the spectral values of Signal B at f = — . 

A secondary indication of the improper behavior of Signal B 
in comparison to atmospheric turbulence can be seen in figure 19. 
Here Signal B was processed by the frequency averaging algorithm 
M = 16 with the Hann window (fig. 19(a)), and with the boxcar 
window (fig. 19(b)). It is apparent that the boxcar window sup- 
pressed random error to a considerably greater extent than did 
the Hann window, an effect exactly opposite from that which would 
be obtained with real turbulence data. A probable explanation is 

— 1 +1 

that the Hann window does not have nulls at — and — as does 

P P 
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the boxcar window. (See fig. 2(b).) The inherent characteristics 
of Signal B would thus appear to make it unsuitable for use in the 
present study. 


Signal A 

Signal A was also generated to contain 16 384 data points 
spaced at an interval of At s 0.04 sec, which yielded a Nyquist 
frequency 12.5 Hz. The signal was constructed as follows: The 

total frequency range, 0 to 12.5 Hz, was partitioned into 1250 
equal intervals of 0.01 Hz each. For each 0.01 Hz interval, a 
uniform-distribution random number generator was used to select 
a frequency f and a phase angle e , with f lying in the inter- 
val and -IT i e ^ IT. For each frequency and phase angle thus 
selected, the sinusoid H(f) sin (2nfn At + e ) (where n = 1, 2, 

. . ., 16 384) was computed. Signal A is the sum of these 1250 

sinusoids. The transfer function H(f) used to adjust the ampli- 
tude of the individual sinusoids was, of course, the same one used 
for Signal B, that is, the square root of equation (6), with 
L = 191 m (627 ft) and V = 183 m/s (600 ft/sec). 

Since both the frequency and phase are made random by the con- 
struction procedure. Signal A is not made up exclusively of har- 
monic components which coincide with boxcar window nulls. The 
approximation of atmospheric turbulence is thus improved. As 
described in the main body of this report, a short segment of Sig- 
nal A is shown in figure 8. The cumulative distribution function, 
shown in figure 9, is very nearly Gaussian, as desired, and as 
expected from the procedure used. 

An additional point worth noting is that the frequency spac- 
ing of the sinusoids, one in each 0.01 Hz interval, is such that 
the smallest bandwidth B employed to process Signal A (about 
0.024 Hz) encompassed at least two sinusoids. Although not tested 
experimentally, the random power fluctuations obtained in the spec- 
tral analysis (for a given B ) could in all likelihood be reduced 
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to an arbitrarily low level by increasing the number of sinusoids 
in each frequency interval. The effect would be similar to that 
achieved by increasing the sample length of actual turbulence data. 
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Figure 3>- Illustration of bias error resulting from convolution of spectral window 

with true spectrum. 
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Figure 7.- Frequency averaging window for M n 16 (sum of 16 squared Hann windows), 
showing bandwidth and window overlap for two adjacent spectral estimates. M, 
number of contiguous raw power spectral density estimates used to compute each 
final estimate at i Af; and K, total number of data points transformed. 
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Figure 8.- Sample of artificially generated time history. Signal A 
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Figure 9.- Cumulative probability distribution of artificially generated 

time history. Signal A. 
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Figure 11.- Power spectrum obtained from ensemble averaging algorithm 
whitening; ^f = 0.0244; B = 0.0366; and B = 0.240. 
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Figure 16.- First-difference prewhitened Dryden spectrum. 
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Figure 17.- Comparison of prewhitened results and non-prewhitened results for wideband 
Blackman-Tukey algorithm, i = 65; Af = 0.195; Bg = 0.391; and ^ei,t " 2.56. 
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